#include "fortran.def"
#include "phys_const.def"
      subroutine cie_thin_cooling_rate(temp, value)
c     written by Tom Abel    8/2003
c     compute optically thin cooling rate due to CIE cooling
c     as discussed in Ripamonti and Abel 2003.
c     input: temp is temperature in Kelvin
c     output total CIE cooling rate (including H2-H2 and He-H2)
C     in erg per second times cm^3 per gram per (gram in H2 molecules).
C     Hence to get cooling rate in erg/s cm^-3 one has to multiply by
C     the total mass density and the mass density in H2 molecules
C     (*rho*rho_H2).
C     This again is an approximation in that it assumes large H2 fractions 
C     (greater than 0.5).
C     Currently it is not clear what happens in the regime where H2
C     becomes mostly dissociated.
      implicit none
#include "fortran_types.def"
      REAL*8 temp
      REAL*8 value
      REAL*8 cie_table(288)
      REAL*8 t_cie(288)
      data t_cie/4.000d+02,4.040d+02,4.080d+02,4.121d+02,4.162d+02,
     $     4.204d+02,4.246d+02,4.289d+02,4.331d+02,4.375d+02,4.418d+02,
     $     4.463d+02,4.507d+02,4.552d+02,4.598d+02,4.644d+02,4.690d+02,
     $     4.737d+02,4.785d+02,4.832d+02,4.881d+02,4.930d+02,4.979d+02,
     $     5.029d+02,5.079d+02,5.130d+02,5.181d+02,5.233d+02,5.285d+02,
     $     5.338d+02,5.391d+02,5.445d+02,5.500d+02,5.555d+02,5.610d+02,
     $     5.666d+02,5.723d+02,5.780d+02,5.838d+02,5.896d+02,5.955d+02,
     $     6.015d+02,6.075d+02,6.136d+02,6.197d+02,6.259d+02,6.322d+02,
     $     6.385d+02,6.449d+02,6.513d+02,6.579d+02,6.644d+02,6.711d+02,
     $     6.778d+02,6.846d+02,6.914d+02,6.983d+02,7.053d+02,7.124d+02,
     $     7.195d+02,7.267d+02,7.339d+02,7.413d+02,7.487d+02,7.562d+02,
     $     7.637d+02,7.714d+02,7.791d+02,7.869d+02,7.948d+02,8.027d+02,
     $     8.107d+02,8.188d+02,8.270d+02,8.353d+02,8.437d+02,8.521d+02,
     $     8.606d+02,8.692d+02,8.779d+02,8.867d+02,8.956d+02,9.045d+02,
     $     9.136d+02,9.227d+02,9.319d+02,9.412d+02,9.506d+02,9.602d+02,
     $     9.698d+02,9.795d+02,9.892d+02,9.991d+02,1.009d+03,1.019d+03,
     $     1.029d+03,1.040d+03,1.050d+03,1.061d+03,1.071d+03,1.082d+03,
     $     1.093d+03,1.104d+03,1.115d+03,1.126d+03,1.137d+03,1.148d+03,
     $     1.160d+03,1.172d+03,1.183d+03,1.195d+03,1.207d+03,1.219d+03,
     $     1.231d+03,1.244d+03,1.256d+03,1.269d+03,1.281d+03,1.294d+03,
     $     1.307d+03,1.320d+03,1.333d+03,1.347d+03,1.360d+03,1.374d+03,
     $     1.387d+03,1.401d+03,1.415d+03,1.430d+03,1.444d+03,1.458d+03,
     $     1.473d+03,1.488d+03,1.502d+03,1.517d+03,1.533d+03,1.548d+03,
     $     1.563d+03,1.579d+03,1.595d+03,1.611d+03,1.627d+03,1.643d+03,
     $     1.660d+03,1.676d+03,1.693d+03,1.710d+03,1.727d+03,1.744d+03,
     $     1.762d+03,1.779d+03,1.797d+03,1.815d+03,1.833d+03,1.852d+03,
     $     1.870d+03,1.889d+03,1.908d+03,1.927d+03,1.946d+03,1.966d+03,
     $     1.985d+03,2.005d+03,2.025d+03,2.045d+03,2.066d+03,2.086d+03,
     $     2.107d+03,2.128d+03,2.150d+03,2.171d+03,2.193d+03,2.215d+03,
     $     2.237d+03,2.259d+03,2.282d+03,2.305d+03,2.328d+03,2.351d+03,
     $     2.375d+03,2.398d+03,2.422d+03,2.447d+03,2.471d+03,2.496d+03,
     $     2.521d+03,2.546d+03,2.571d+03,2.597d+03,2.623d+03,2.649d+03,
     $     2.676d+03,2.702d+03,2.730d+03,2.757d+03,2.784d+03,2.812d+03,
     $     2.840d+03,2.869d+03,2.897d+03,2.926d+03,2.956d+03,2.985d+03,
     $     3.015d+03,3.045d+03,3.076d+03,3.106d+03,3.138d+03,3.169d+03,
     $     3.201d+03,3.233d+03,3.265d+03,3.298d+03,3.331d+03,3.364d+03,
     $     3.397d+03,3.431d+03,3.466d+03,3.500d+03,3.535d+03,3.571d+03,
     $     3.606d+03,3.643d+03,3.679d+03,3.716d+03,3.753d+03,3.790d+03,
     $     3.828d+03,3.867d+03,3.905d+03,3.944d+03,3.984d+03,4.024d+03,
     $     4.064d+03,4.105d+03,4.146d+03,4.187d+03,4.229d+03,4.271d+03,
     $     4.314d+03,4.357d+03,4.401d+03,4.445d+03,4.489d+03,4.534d+03,
     $     4.579d+03,4.625d+03,4.671d+03,4.718d+03,4.765d+03,4.813d+03,
     $     4.861d+03,4.910d+03,4.959d+03,5.008d+03,5.058d+03,5.109d+03,
     $     5.160d+03,5.212d+03,5.264d+03,5.316d+03,5.370d+03,5.423d+03,
     $     5.477d+03,5.532d+03,5.588d+03,5.643d+03,5.700d+03,5.757d+03,
     $     5.814d+03,5.873d+03,5.931d+03,5.991d+03,6.051d+03,6.111d+03,
     $     6.172d+03,6.234d+03,6.296d+03,6.359d+03,6.423d+03,6.487d+03,
     $     6.552d+03,6.617d+03,6.684d+03,6.750d+03,6.818d+03,6.886d+03,
     $     6.955d+03/
      data cie_table/1.435d-15,1.488d-15,1.544d-15,1.601d-15,1.661d-15,
     $     1.722d-15,1.786d-15,1.853d-15,1.922d-15,1.993d-15,2.067d-15,
     $     2.143d-15,2.223d-15,2.305d-15,2.390d-15,2.479d-15,2.570d-15,
     $     2.665d-15,2.763d-15,2.865d-15,2.971d-15,3.080d-15,3.194d-15,
     $     3.314d-15,3.440d-15,3.571d-15,3.707d-15,3.848d-15,3.995d-15,
     $     4.148d-15,4.307d-15,4.472d-15,4.643d-15,4.822d-15,5.008d-15,
     $     5.201d-15,5.402d-15,5.611d-15,5.829d-15,6.056d-15,6.292d-15,
     $     6.539d-15,6.800d-15,7.071d-15,7.354d-15,7.648d-15,7.955d-15,
     $     8.274d-15,8.607d-15,8.954d-15,9.316d-15,9.692d-15,1.008d-14,
     $     1.049d-14,1.092d-14,1.136d-14,1.183d-14,1.235d-14,1.291d-14,
     $     1.350d-14,1.412d-14,1.477d-14,1.545d-14,1.617d-14,1.693d-14,
     $     1.773d-14,1.857d-14,1.946d-14,2.039d-14,2.137d-14,2.242d-14,
     $     2.352d-14,2.469d-14,2.592d-14,2.721d-14,2.857d-14,3.000d-14,
     $     3.151d-14,3.310d-14,3.477d-14,3.653d-14,3.838d-14,4.053d-14,
     $     4.300d-14,4.561d-14,4.837d-14,5.127d-14,5.434d-14,5.758d-14,
     $     6.099d-14,6.458d-14,6.837d-14,7.236d-14,7.687d-14,8.167d-14,
     $     8.675d-14,9.210d-14,9.776d-14,1.037d-13,1.100d-13,1.167d-13,
     $     1.237d-13,1.310d-13,1.388d-13,1.470d-13,1.556d-13,1.646d-13,
     $     1.741d-13,1.841d-13,1.946d-13,2.057d-13,2.172d-13,2.294d-13,
     $     2.421d-13,2.554d-13,2.694d-13,2.841d-13,2.994d-13,3.155d-13,
     $     3.323d-13,3.498d-13,3.682d-13,3.874d-13,4.074d-13,4.284d-13,
     $     4.502d-13,4.730d-13,4.968d-13,5.216d-13,5.475d-13,5.744d-13,
     $     6.025d-13,6.318d-13,6.622d-13,6.939d-13,7.268d-13,7.611d-13,
     $     7.967d-13,8.337d-13,8.721d-13,9.121d-13,9.535d-13,9.966d-13,
     $     1.041d-12,1.088d-12,1.136d-12,1.185d-12,1.237d-12,1.290d-12,
     $     1.346d-12,1.403d-12,1.462d-12,1.524d-12,1.587d-12,1.653d-12,
     $     1.721d-12,1.791d-12,1.864d-12,1.939d-12,2.016d-12,2.096d-12,
     $     2.179d-12,2.268d-12,2.374d-12,2.483d-12,2.597d-12,2.714d-12,
     $     2.836d-12,2.961d-12,3.091d-12,3.226d-12,3.365d-12,3.508d-12,
     $     3.657d-12,3.810d-12,3.968d-12,4.132d-12,4.300d-12,4.474d-12,
     $     4.654d-12,4.839d-12,5.029d-12,5.226d-12,5.428d-12,5.637d-12,
     $     5.852d-12,6.073d-12,6.301d-12,6.535d-12,6.776d-12,7.024d-12,
     $     7.280d-12,7.542d-12,7.812d-12,8.089d-12,8.374d-12,8.666d-12,
     $     8.967d-12,9.276d-12,9.593d-12,9.919d-12,1.025d-11,1.060d-11,
     $     1.099d-11,1.145d-11,1.193d-11,1.241d-11,1.292d-11,1.343d-11,
     $     1.396d-11,1.451d-11,1.507d-11,1.565d-11,1.624d-11,1.685d-11,
     $     1.748d-11,1.813d-11,1.879d-11,1.947d-11,2.017d-11,2.089d-11,
     $     2.162d-11,2.238d-11,2.316d-11,2.395d-11,2.477d-11,2.561d-11,
     $     2.646d-11,2.734d-11,2.824d-11,2.917d-11,3.011d-11,3.118d-11,
     $     3.234d-11,3.353d-11,3.475d-11,3.601d-11,3.729d-11,3.861d-11,
     $     3.996d-11,4.135d-11,4.277d-11,4.422d-11,4.571d-11,4.724d-11,
     $     4.881d-11,5.041d-11,5.205d-11,5.373d-11,5.545d-11,5.721d-11,
     $     5.901d-11,6.085d-11,6.273d-11,6.470d-11,6.692d-11,6.920d-11,
     $     7.152d-11,7.391d-11,7.635d-11,7.885d-11,8.141d-11,8.402d-11,
     $     8.670d-11,8.944d-11,9.224d-11,9.510d-11,9.803d-11,1.010d-10,
     $     1.041d-10,1.072d-10,1.104d-10,1.137d-10,1.170d-10,1.204d-10,
     $     1.239d-10,1.275d-10,1.311d-10,1.348d-10,1.386d-10,1.425d-10,
     $     1.464d-10,1.505d-10,1.546d-10,1.588d-10,1.631d-10,1.674d-10,
     $     1.719d-10/
      INTG_PREC i, j, jmin, jmax
      REAL*8 mh
      mh = mass_h
c     rough low temperature extrapolation
      if (temp .le. t_cie(1)) then
         value = cie_table(1)*(temp/t_cie(1))**4
c            print*, value
c         write(6,*), "inside_cie: ", temp, value
         return
      endif
c     extrapolate high temperature with **3. Not discussed in Ripamonti
c     & Abel ...
      if (temp .ge. t_cie(288)) then
         value = cie_table(288)*(temp/t_cie(288))**3
c            print*, value
c         write(6,*), "inside_cie: ", temp, value
         return
      endif
      jmin = 1 
      jmax = 288
      do i=0,200
         j = (jmin+jmax)/2
         if (temp .ge. t_cie(j)) then 
            jmin=j
         else 
            jmax=j
         endif            
         if ((jmax-jmin) .le. 1) then
c            write(6,*) j, jmin, jmax,cie_table(jmin),cie_table(jmax)
c            write(6,*) t_cie(jmin),temp,t_cie(jmax)
            value =
     $           (cie_table(jmax)*(temp-t_cie(jmin))+
     $            cie_table(jmin)*(t_cie(jmax)-temp))/
     $            (t_cie(jmax)-t_cie(jmin))
c            write(6,*), "inside_cie: ", temp, value
            return
         endif
      enddo
      return 
      end
